#' -----------------------------------------------
#' 
#' Code for:
#' 
#' Politicians’ misinformation, its correction, 
#' and partisanship in Italy
#' 
#' Authors:
#' Mario Quaranta
#' Luca M. Arrigoni
#' 
#' Published in: 
#' Italian Political Science Review /
#' Rivista Italiana di Scienza Politica
#'
#' -----------------------------------------------

# Libraries -----

library(tidyverse)

# Set wd ----

setwd("/~")

# Load dataset -----

load(file = "df.Rdata")

# Functions ----

extr_preds <- function(model, term){
  
  eff <- effects::effect(term, model)
  
  X <- eff$model.matrix
  
  S <- MASS::mvrnorm(100000, coef(model), vcov(model))
  
  y <- S %*% t(X)
  colnames(y) <- c("control", "ec1", "ec2", "ec3")
  
  ci <- apply(y, 2, quantile, c(0.025, .5, .975))
  colnames(ci) <- c("control", "ec1", "ec2", "ec3")
  
  diff_control <- apply(y[, 1:4] - y[, "control"], 2, quantile, c(0.025, .5, .975))
  diff_ec1 <- apply(y[, 2:4] - y[, "ec1"], 2, quantile, c(0.025, .5, .975))
  
  return(list(ci = ci, diff_control = diff_control, diff_ec1 = diff_ec1))
  
}

extr_preds_int <- function(model, term, levs = levs, xlevels = xlevels){
  
  eff <- effects::effect(term, model, xlevels = xlevels)
  
  X <- eff$model.matrix
  
  S <- MASS::mvrnorm(100000, coef(model), vcov(model))
  
  y <- S %*% t(X)
  
  colnames(y) <- paste(rep(paste("EC", 0:3), 3), "-", rep(levs, each = 4)) 
  
  ci <- apply(y, 2, quantile, c(0.025, .5, .975))
  
  ci <- t(ci)
  ci <- as.data.frame(ci)
  ci$v <- factor(paste("EC", 0:3))
  ci$x <-rep(levs, each = 4)
  
  ci <- ci %>% arrange(v, x)
  names(ci)[1:3] <- c("low", "y", "upp")
  
  start <- seq(1, ncol(y), by = 4)
  
  sdf <- lapply(start, function(i, y) y[, i:(i + 3)], y = y)
  
  diff <- sdf
  
  for(i in 1:length(levs)){
    diff[[i]] <- sdf[[i]][, 1:4] - sdf[[i]][, 1]
  }
  
  diff <- lapply(diff, function(x) t(apply(x, 2, quantile, c(0.025, .5, .975))))
  diff <- as.data.frame(do.call(rbind, diff))
  diff$x <- rep(levs, each = 4)
  diff$c <- 0:3
  
  return(list(ci = ci, eff = eff, diff = diff))
  
}

# Descriptive statistics (for appendix) ----

df_desc_cont <- df %>% 
  select(
    dv_exp1, dv_exp2, dv_exp3,
    feel_pd, feel_fdi, feel_m5s,
    Q46
  )

df_desc_cat <- df %>% 
  select(
    ballot_1, ballot_2, ballot_3, 
    Q45, Q47
  )

cont_m <- apply(df_desc_cont, 2, mean, na.rm = T)
cont_s <- apply(df_desc_cont, 2, sd, na.rm = T)
cont_n <- apply(df_desc_cont, 2, function(x) sum(!is.na(x)))

props <- apply(df_desc_cat, 2, function(x) as.data.frame(prop.table(table(x))))
freq  <- apply(df_desc_cat, 2, function(x) as.data.frame(table(x)))

contain <- rep(NA, ncol(df_desc_cat))

for(i in 1:ncol(df_desc_cat)){
  contain[i] <- sum(!is.na(df_desc_cat[, i]))
}

props <- do.call(rbind, props)
freq  <- do.call(rbind, freq)

cont <- cbind(cont_m, cont_s, cont_n)
colnames(cont) <- c("m", "s", "n")

prop <- cbind(props[, 2], rep(NA, nrow(freq)), freq[, 2])
colnames(prop) <- c("m", "s", "n")
rownames(prop) <- props[, 1]

desc_tab <- rbind(cont, prop)

desc_tab <- desc_tab[c(1:3, 8:19, 4:6, 20:21, 7, 23, 24, 22), ]

desc_tab[, c(1, 2)] <- format(round(desc_tab[, c(1, 2)], 3), 4)

write.table(desc_tab, file = "tab_b1.txt", sep = "\t")

# Models (baseline, for appendix) ----

fit_1_bl <- lm(dv_exp1 ~ ballot_1, data = df)
fit_2_bl <- lm(dv_exp2 ~ ballot_2, data = df)
fit_3_bl <- lm(dv_exp3 ~ ballot_3, data = df)

# Models (w/ controls) ----

fit_1 <- lm(dv_exp1 ~ ballot_1 + feel_pd + Q45 + Q46 + Q47, data = df)
fit_2 <- lm(dv_exp2 ~ ballot_2 + feel_fdi + Q45 + Q46 + Q47, data = df)
fit_3 <- lm(dv_exp3 ~ ballot_3 + feel_m5s + Q45 + Q46 + Q47, data = df)

# Extract values for plots

y_1 <- extr_preds(fit_1, "ballot_1")
y_2 <- extr_preds(fit_2, "ballot_2")
y_3 <- extr_preds(fit_3, "ballot_3")

# Plot predicted values (for appendix)

df_plot <- cbind(y_1$ci, 
                 y_2$ci,
                 y_3$ci)

df_plot <- t(df_plot)
df_plot <- as.data.frame(df_plot)
names(df_plot) <- c("low", "y", "upp")
df_plot$politician <- rep(1:3, each = 4)
df_plot$politician <- factor(df_plot$politician, 1:3, 
                             c("PD / Minimum wage", "FDI / Vaccination campaign", "M5S / Working hours"))

df_plot$condition <- 1:4
df_plot$condition <- factor(df_plot$condition, 1:4, c("Control", paste("EC", 1:3)))

df_plot %>% 
  ggplot(aes(x = condition, y = y)) + 
  geom_point(position = position_dodge(width = 0.3)) + 
  geom_segment(aes(y = low, yend = upp), position = position_dodge(width = 0.3)) +
  
  facet_grid(~ politician) +
  labs(x = "", y = "Predicted values") +
  scale_y_continuous(limits = c(3, 9), breaks = seq(2, 10, 1), labels = seq(2, 10, 1)) + 
  theme_light() +
  theme(strip.background = element_rect(fill = "white"),
        strip.text = element_text(color = "black"),
        axis.text.x = element_text(size = 7))


ggsave(file = "figure-a1.pdf", width = 7, height = 3)
ggsave(file = "figure-a1.png", width = 7, height = 3, dpi = 1200)

# Plot all groups vs control (for article)

df_plot <- cbind(y_1$diff_control, 
                 y_2$diff_control,
                 y_3$diff_control 
)

df_plot <- t(df_plot)
df_plot <- as.data.frame(df_plot)
names(df_plot) <- c("low", "y", "upp")
df_plot$politician <- rep(1:3, each = 4)
df_plot$politician <- factor(df_plot$politician, 1:3, 
                             c("PD / Minimum wage", "FDI / Vaccination campaign", "M5S / Working hours"))

df_plot$condition <- 1:4
df_plot$condition <- factor(df_plot$condition, 1:4, c("Control", paste("EC", 1:3)))

df_plot %>% 
  filter(condition != "Control") %>% 
  ggplot(aes(x = condition, y = y)) + 
  geom_point(position = position_dodge(width = 0.3)) + 
  geom_segment(aes(y = low, yend = upp), position = position_dodge(width = 0.3)) +
  geom_hline(yintercept = 0, linetype = 2) + 
  facet_grid(~ politician) +
  labs(x = "", y = "Predicted differences") +
  scale_y_continuous(limits = c(-1.5, 1.5), breaks = seq(-1.5, 1.5, .5), labels = seq(-1.5, 1.5, .5)) + 
  theme_light() +
  theme(strip.background = element_rect(fill = "white"),
        strip.text = element_text(color = "black"),
        axis.text.x = element_text(size = 7))

ggsave(file = "figure-1.pdf", width = 7, height = 3)
ggsave(file = "figure-1.png", width = 7, height = 3, dpi = 1200)

# Differences EC2 vs. EC1 - EC3 vs. EC1 (for appendix)

df_plot <- cbind(y_1$diff_ec1, 
                 y_2$diff_ec1,
                 y_3$diff_ec1
)

df_plot <- as.data.frame(t(df_plot))
names(df_plot) <- c("low", "y", "upp")
df_plot$condition <- 1:3
df_plot$condition <- factor(df_plot$condition, 1:3, c("EC 1", "EC 2", "EC 3"))
df_plot$politician <- rep(1:3, each = 3)
df_plot$politician <- factor(df_plot$politician, 1:3, 
                             c("PD / Minimum wage", "FDI / Vaccination campaign", "M5S / Working hours"))

df_plot %>% 
  filter(condition != "EC 1") %>% 
  ggplot(aes(x = condition, y = y)) + 
  geom_point(position = position_dodge(width = 0.3)) + 
  geom_segment(aes(y = low, yend = upp), position = position_dodge(width = 0.3)) +
  
  geom_hline(yintercept = 0, linetype = 2) + 
  facet_grid(~ politician) +
  labs(x = "", y = "Predicted differences") +
  scale_y_continuous(limits = c(-1.5, 1.5), breaks = seq(-1.5, 1.5, .5), labels = seq(-1.5, 1.5, .5)) + 
  theme_light() +
  theme(strip.background = element_rect(fill = "white"),
        strip.text = element_text(color = "black"),
        axis.text.x = element_text(size = 7))

ggsave(file = "figure-a2.pdf", width = 7, height = 3)
ggsave(file = "figure-a2.png", width = 7, height = 3, dpi = 1200)

# Models table (for appendix)

tab_est <- cbind(
  rbind(summary(fit_1_bl)$coefficients[, c(1, 2, 4)],
        matrix(NA, ncol = 3, nrow = 5)),  
  summary(fit_1)$coefficients[, c(1, 2, 4)],
  rbind(summary(fit_2_bl)$coefficients[, c(1, 2, 4)],
        matrix(NA, ncol = 3, nrow = 5)),  
  summary(fit_2)$coefficients[, c(1, 2, 4)],
  rbind(summary(fit_3_bl)$coefficients[, c(1, 2, 4)],
        matrix(NA, ncol = 3, nrow = 5)),  
  summary(fit_3)$coefficients[, c(1, 2, 4)]
)

tab_est <- format(round(tab_est, 3), digits = 4)

tab_est <- rbind(tab_est[1, ],
                 rep(NA, 18),
                 tab_est[2:7, ],
                 rep(NA, 18),
                 tab_est[8:9, ],
                 c(length(fit_1_bl$fitted.values), NA, NA, 
                   length(fit_1$fitted.values), NA, NA,
                   length(fit_2_bl$fitted.values), NA, NA,
                   length(fit_2$fitted.values), NA, NA,
                   length(fit_3_bl$fitted.values), NA, NA,
                   length(fit_3$fitted.values), NA, NA))

rownames(tab_est) <- c("Intercept", "Group (r.c. Control):", paste("EC", 1:3),
                       "Feeling t/w party", "Gender (Man)", 
                       "Age", "Education (r.c. Low):", "Medium", "High",
                       "N")

tab_b2 <- tab_est[c(1:5, 12), c(1:3, 7:9, 13:15)]
tab_b3 <- tab_est[, c(4:6, 10:12, 16:18)]

write.table(tab_b2, file = "tab_b2.txt", sep = "\t")
write.table(tab_b3, file = "tab_b3.txt", sep = "\t")

# Models (w/ interactions) ----

fit_1_int <- lm(dv_exp1 ~ ballot_1 * feel_pd +  Q45 + Q46 + Q47, data = df)
fit_2_int <- lm(dv_exp2 ~ ballot_2 * feel_fdi + Q45 + Q46 + Q47, data = df)
fit_3_int <- lm(dv_exp3 ~ ballot_3 * feel_m5s + Q45 + Q46 + Q47, data = df)

# Extract values for plots

levs <- 0:10
y_1_int <- extr_preds_int(fit_1_int, "ballot_1:feel_pd",  levs = levs, xlevels = list(feel_pd = levs))
y_2_int <- extr_preds_int(fit_2_int, "ballot_2:feel_fdi", levs = levs, xlevels = list(feel_fdi = levs))
y_3_int <- extr_preds_int(fit_3_int, "ballot_3:feel_m5s", levs = levs, xlevels = list(feel_m5s = levs))

# Predicted values (for appendix)

df_plot <- rbind(y_1_int$ci, 
                 y_2_int$ci,
                 y_3_int$ci
)

names(df_plot)[1:3] <- c("low", "y", "upp")
df_plot$politician <- rep(1:3, each = nrow(y_1_int$ci))
df_plot$politician <- factor(df_plot$politician, 1:3, 
                             c("PD / Minimum wage", "FDI / Vaccination campaign", "M5S / Working hours"))

df_plot$x <- factor(df_plot$x)

df_plot$c <- factor(df_plot$v, paste("EC", 0:3), c("Control", paste("EC", 1:3)))

df_plot %>% 
  filter(x %in% c("0", "5", "10")) %>% 
  ggplot(aes(x = c, y = y, color = x, shape = x)) + 
  geom_segment(aes(y = low, yend = upp), position = position_dodge(width = 0.3)) +
  geom_point(position = position_dodge(width = 0.3)) +
  facet_grid( ~ politician) +
  labs(x = "", y = "Predicted values", color = "Feeling", shape = "Feeling") +
  scale_color_manual(values = c("black", "grey40", "grey")) +
  scale_y_continuous(limits = c(3, 9), breaks = seq(2, 10, 1), labels = seq(2, 10, 1)) + 
  theme_light() +
  theme(
    strip.background = element_rect(fill = "white"),
    strip.text = element_text(color = "black"),
    legend.title = element_text(size = 8),
    axis.text.x = element_text(size = 8),
    legend.position = "top") + guides(fill = "none")

ggsave(file = "figure-a3.pdf", width = 8.5, height = 4)
ggsave(file = "figure-a3.png", width = 8.5, height = 4, dpi = 1200)

# Plot all groups vs control (for article)

df_plot <- rbind(y_1_int$diff, 
                 y_2_int$diff,
                 y_3_int$diff
)

names(df_plot)[1:3] <- c("low", "y", "upp")
df_plot$politician <- rep(1:3, each = nrow(y_1_int$diff))
df_plot$politician <- factor(df_plot$politician, 1:3, 
                             c("PD / Minimum wage", "FDI / Vaccination campaign", "M5S / Working hours"))

df_plot$x <- factor(df_plot$x)

df_plot$c <- factor(df_plot$c, 0:3, c("Control", paste("EC", 1:3)))

df_plot %>% 
  filter(x %in% c("0", "5", "10"),
         c != "Control") %>% 
  ggplot(aes(x = c, y = y, color = x, shape = x)) + 
  geom_segment(aes(y = low, yend = upp), position = position_dodge(width = 0.3)) +
  geom_point(position = position_dodge(width = 0.3)) +
  facet_grid( ~ politician) +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(x = "", y = "Predicted differences", color = "Feeling", shape = "Feeling") +
  scale_color_manual(values = c("black", "grey40", "grey")) +
  scale_y_continuous(limits = c(-2.0, 2.6), breaks = seq(-2, 2.5, .5), labels = seq(-2, 2.5, .5)) + 
  theme_light() +
  theme(strip.background = element_rect(fill = "white"),
        strip.text = element_text(color = "black"),
        legend.title = element_text(size = 8),
        axis.text.x = element_text(size = 8),
        legend.position = "top") + 
  guides(fill = "none")

ggsave(file = "figure-2.pdf", width = 8.5, height = 4)
ggsave(file = "figure-2.png", width = 8.5, height = 4, dpi = 1200)

# Models table (for appendix)

tab_est <- cbind(
  summary(fit_1_int)$coefficients[, c(1:2, 4)],
  summary(fit_2_int)$coefficients[, c(1:2, 4)],
  summary(fit_3_int)$coefficients[, c(1:2, 4)]
)

tab_est <- format(round(tab_est, 3), digits = 4)

tab_est <- rbind(tab_est[1, ],
                 rep(NA, 9),
                 tab_est[2:7, ],
                 rep(NA, 9),
                 tab_est[8:9, ],
                 rep(NA, 9),
                 tab_est[10:12, ],
                 c(length(fit_1_int$fitted.values), NA, NA,
                   length(fit_2_int$fitted.values), NA, NA,
                   length(fit_3_int$fitted.values), NA, NA))

rownames(tab_est) <- c("Intercept", "Group (r.c. Control):", paste("EC", 1:3),
                       "Feeling t/w party", "Gender (Man)", 
                       "Age", "Education (r.c. Low):", "Medium", "High",
                       "Interaction",
                       paste(paste("EC", 1:3), "* Feeling t/w party"),
                       "N")

write.table(tab_est, file = "tab_b4.txt", sep = "\t")
